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ABSTRACT 

We study the gravitational instability (GI) of small solids in a gas disk as a mechanism to form 
planetesimals. Dissipation from gas drag introduces secular GI, which proceeds even when standard 
GI criteria for a critical density or Toomre's Q predict stability. We include the stabilizing effects of 
turbulent diffusion, which suppresses small scale GI. The radially wide rings that do collapse contain 
up to ^ 0.1 Earth masses of solids. Subsequent fragmentation of the ring (not modeled here) would 
produce a clan of chemically homogenous planetesimals. Particle radial drift time scales (and, to a 
lesser extent, disk lifetimes and sizes) restrict the viability of secular GI to disks with weak turbulent 
diffusion, characterized by a < 10""*. Thus midplane dead zones are a preferred environment. Large 
solids with radii > 10 cm collapse most rapidly because they partially decouple from the gas disk. 
Smaller solids, even below ^ mm-sizes could collapse if particle-driven turbulence is weakened by 
either localized pressure maxima or super-Solar metallicity. Comparison with simulations that include 
particle clumping by the streaming instability shows that our linear model underpredicts rapid, small 
scale gravitational collapse. Thus the inclusion of more detailed gas dynamics promotes the formation 
of planetesimals. We discuss relevant constraints from Solar System and accretion disk observations. 
Subject headings: hydrodynamics — instabilities — planetary systems: formation — planetary sys- 
tems: protoplanetary disks — solar system: formation — turbulence 



1. INTRODUCTION 

The discovery of a diverse population of extrasolar 
planet^ motivates the development of planet formation 
theories that are (at least) physically consistent and (at 
best) predictive of the rapidly growing data set. Form- 
ing planetesimals above kilometer sizes is both a cru- 
cial first step an d testa ble in our own Solar system. 
[Chiang & Y^dii pOlol hereafter CYIO) and lYoudinl 
^010) review the main physical processes and theoretical 
puzzles. 

A promising route to planetesimal formation is 
the clumping of solids by the streaming instabil- 
ity ( SI ) in protoplanetary disks (lYoudin &: Goodnian 
2005 : lYoudin fc JohansenI 120071: iJohansen fc Youdin 
1007), wh ich then induces gravitational collapse 
Johansen _ et al. 2007 . 2009b'). SI is a drag instability 



Goodman fc Pindoill2000i ) that arises from the aerody- 



namic back reaction of solids on gas, and does not require 
self-gravity. 

Strong Sl-induced clumping appears to require fairly 
large solids with Ts ~ 1, where 

Ts = /^tstop (1) 

measures the dynamical relevance of drag forces with 
tstop the time scale for damping particle relative mo- 
tion and f2 the Keplerian orbital frequency. Thus Tg 
increases with particle radius a. In standard disk mod- 
els, marginal coupling, Ts = 1, corresponds to a < 1 
meter near 1 AU, and smaller sizes in the outer disk 
where gas densities are lower (see Fig. [T|). For small 
solids with Tg ^ 1, SI operates but it stirs particles 



without producing strong clumping (jJohansen fc Youdlnl 
I2007t IBai fc Stonell2010aj) . Most other particle concen- 
tration me chanisms also become inefficient for Tg <C 1 
(see CYlO. r^udh3[2OT0l) . 

It is desirable to find a planetesimal formation mech- 
anism that operates when Ts <^ 1. Protoplanetary 
disks initially contain sub-micron-sized grains. Col- 
lisional sticking alo ne ap pears un likely to p r oduce 
meter-sized bodies (iBma 12 000: Y oudin fc Shul [20021: 
ISekiva fc Takedal 120031: 1 Yoiodini, 2004: IZsom et al.l[2oTol) . 
Furthermore, many primitive meteorites are composed 
prim arily of chondrules, in clusions with a ^ 0.1 — 1 
mm ([Russeh et all [2001 . iCuzzi et all pOOl propose 
that turbulent eddies concentrate chondrules into clumps 
that become 10 — 100 km planetesimals. As discussed 
more extensively in CYIO, this concentration mechanism 
normally applies to much smaller scales. More study of 
mechanisms to concentrate Tg ^ 1 solids is needed. 

This work considers gravitational instabilities (GI) 
that are secular, meaning mediated by dissipation in the 
form of gas drag. T he original proposals that planetes- 
imals formed by GI (ISafronovl 119691: IGoldreich fc WaTdl 



19731) applied t he standard, dissipationless GI criterion of 



J iiiea t r 
1964f ) . This approach gives the classic result of 



I ayoudin@cf'a. harvard, edu | 

^ See http://exoplanet.eu', compiled by Jean Schneider, for an 
updated catalog with reterences. 



~ kilometer-sized planetesimals. The fact that gas drag 
affects GI should not be surprising, especially for Tg ^ 1. 
Section 16.51 summarizes previous work on secular GI . 

This paper contains an improved treatment of gas 
turbulence, notably the stabilizing effects of turbu- 
lent diffusion. We use th e turbulent stirring model of 
lYoudin fc LithwicH (|2007[ ). which includes the effects of 
orbital oscillations. Our simplified treatment of gas dy- 
namics does not capture SI or other non-gravitational 
concentration mechanisms. However, we do vary the disk 
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"metallicity" 



(2) 



the ratio of particle (unsubscripted) to gas surface densi- 
ties. Large Z values mimic particle concentration or gas 
depletion on large scales (compared to the wavelength) 
and accelerate secular GI. 

We emphasize one general result from our analysis. 
The critical or Roche-lik£0 density 



(3) 



is not a useful discriminant for gravitational collapse 
when drag forces are significant, as they are for planetes- 
imal formation. We use M* for the central stellar mass 
and R for the radial distance from that star. Dissipation 
allows collapse to proceed for particle densities p <^ p^i- 
We prove this result in !j3.1.3l and further demonstrate 
it in Fig. [6l Moreover, \i p > pR in only a small vol- 
ume, then collapse is avoidable. Gas drag limits collapse 
speeds to the terminal velocity, and turbulent diffusion 
smoothes small scale density perturbations. 

Since secular GI can require many orbits, it is re- 
stricted by radial drift speeds and ultimately disk life- 
times. Observations find protoplanetary disks wi th ages 
up to several million years ([Hernandez et al.l2007| ) . Plan- 
etesimal formation within 1 Myr should still allow time 
for the growth of gas and i ce giants by core accretion 
(jLissauer &: StevensonI I2007D . The existence of prim- 
itive meteoritic inclusion s with age spreads of several 
Myr ([Russell et al.l 120061 ) supports a longer time scale 
for planetesimal formation in the inner Solar System. 

This article is organized as follows. We start with an 
order of magnitude derivation of secular GI in We 
present our full model in f|3] including: a derivation of the 
dispersion relation in i i3.1l the turbulent stirring model in 
§3.21 and the protoplanetary disk model in i i3.3l Section 
U gives the main numerical results covering: the varia- 
tion of GI properties with disk parameters in i j4.2l and 
a calculation of the strongest levels of turbulence that 
allow collapse in ij4.3l Section [5] explores the analytic 
behavior of secular secular GI and the transition to stan- 
dard GI. We discuss several aspects of o ur re sults in ^j6] 
including a comparison to simulations in ^ 36.11 Solar Sys- 
tem constraints in ! j6.2l and metallicity thresholds in !j6.3l 
Section [7] concludes with a brief summary. Appendix |^ 
addresses turbulent drag laws and the behavior of secular 
GI in less relevant regions of parameter space is relegated 
to Appendix IbI 

2. SECULAR COLLAPSE OF SMALL SOLIDS 

We begin with an order-of-magnitude derivation of the 
main result of this paper: the gravitational collapse of 
small solids in a gas disk. Small means that drag forces 
are strong with Tg ^ 1 . Though probably less relevant for 
planetesimal formation, secular GI also exists for Tg 3> 1, 
see gEl 

Consider a ring-like surface density perturbation of am- 
plitude CT and radial length scale A. The radial gravity 

^ The actual Roche density refers to tidal disruption of a sphere 
and is numerically larger (see CYIO). Th e numerical co efficient in 
equation (jsjl is fro m the buckling modes of lSekival l|1983l ). which we 
discuss further in i]6.5l 
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Figure 1. Aerodynamic stopping times (normalized by the or- 
bital frequency) vs. radial distance in the disk midplane. For our 
reference minimum mass model {F = 1), we plot Ts for particle 
sizes: a = 1 mm {black solid curve), a = 3 cm {green dot-dashed 
curve) and a = 30cm, {red dotted curve). For a gas depleted disk 
{F = 0.1) we again plot Ts for a = 3 cm {purple tripe dot-dashed 
curve) and a = 30 cm {blue dashed curve). The kinks in the c urves 
correspond to transitions between drag laws, as described in i|3.3.2l 
and appendix lAl 



at the edge of the perturbation is gn ~ Ga. If this re- 
sult is unfamiliar, recall that the acceleration towards an 
infinitely long wire of linear mass density A is 2GA/X at 
a distance A. Then identify cr ^ A/A as if the ring were 
concentrated in a thin wire, which can be approximated 
as straight and infinitely long because A <C i?, the radial 
distance to the star. 

The particles' response to self-gravity is mediated by 
gas drag. Radial sedimentation to the center of the ring 
occurs at the terminal velocity 

" = gntstop ~ Gatstop ■ (4) 

It is reasonable to ask why the radial gravity induces ra- 
dial collapse instead of merely shifting the orbital speed 
by adjusting centripetal balance. The reason is that 
the drag on this extra orbital motion would would give 
radial infall that exceeds the terminal velocity limit, 
hence equation (jlj applies. This explanation assumes the 
pressure-supported gas remains stationary in response to 
the particle perturbations, a critical issue that we discuss 
further in ^6.5\ 

By mass continuity, the amplitude of the density per- 
turbation grows at a rate 



Su 



GSt 



stop 



A 



(5) 



To understand whether and how fast collapse occurs, we 
must consider stabilizing infiuences that restrict the al- 
lowable A. The excess angular momentum in a pertur- 
bation is removed very rapidly, in igtop ^ 1/s as we can 
verify shortly. Thus there is no obstacle to long wave- 
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length collapse, and instability is guaranteed. We now 
consider impediments to small scale collapse that set the 
minimum A. 

2.1. Effective Pressure Support: Less Relevant 

We first consider the stabilization by the effective pres- 
sure of particle random velocities, c. Gravity, gn exceeds 
the effective pressure acceleration, c'^a/{X!\), for wave- 
lengths above 



c 



(6) 



From equation ([S]) the fastest growth rate allowed by 
effective pressure support is 



sp 



stop ■ 



[2- 



(7) 



Qt = cn/{TTGE) is the standard iToomrd (|1964D gravi- 
tational stability parameter, but unlike standard GI we 
have no requirement that Qt < 1- 

2.2. Diffusive Support: Dominant 

Next we consider the stabilizing effect of particle diffu- 
sivity, Z?, due to stirring by gas turbulence. We express 



D : 



(8) 



as the product of the particle speed, c, and path length, 
c//2, which arises because disk turbulence randomizes 
the particle trajectory in an eddy turnover time toddy — 
1/ n. See §3.21 for more details on turbulent stirring. 

For the collapse rate of equation ([5]) to outpace the 
diffusive spreading time scale A^/I?, wavelengths must 
exceed 

^--ci — 7^- 

Thus \d ^ Ap, by our assumption of strong drag. Dif- 
fusion is the more effective stabilizing influence that sets 
the actual minimum A. 
Secular GI thus has a collapse rate 



(G17tstop)^ „ 



(10) 



that is slower than sp by Tg ^ 1. Growth becomes sig- 
niflcantly slower as particle size, and thus Ts, declines. 

Previous studies o f dissipative collapse that neglect ed 
turbulent diffusion (|Wardlll97fil ^2Mk lYoudii] 120051 0) 
were thus overly optimistic. However the qualitative re- 
sult that collapse does not require Qt < 1 remains. 

Armed with this basic description of the instability, 
some readers may wish to skip ahead to the numerical 
results of 21 and the figures. To better understand the 
model ingredient and caveats, we do recommend coming 
back to ^131 

3. FORMULATION OF THE PROBLEM 

3.1. Dynamical Equations 

We consider the local, linear dynamics of solid particles 
in a protoplanetary disk subject to gas drag, turbulent 
mixing and self-gravity. We adopt the local shearing box 
model for a patch of the disk (,Goldreich k, Lvndcn-Beiil 



Il965f) with a frame centered on a fiducial radius that 
rotates with the local Kepler frequency Q. The lo- 
cal radial and azimuthal coordinates are x and y, re- 
spectively. We decompose the particle surface density 
Ep = E -|- cr and velocity fields V = Vq -\- ux + vy into 
time-steady backgrounds — constant E and the Keple- 
rian shear Vq = — (3/2)J7xy — plus small amplitude 
perturbations. The perturbations are axisymmetric as 
appropriate for slow motions stretched by shear. Evolu- 
tion in time, i, follows the linearized, height-integrated 
equations of continuity and force: 



da ^du 
dt dx 
du 
'dt 

dv ^ ,„ 

- + nui2 



= D 



2Qv 



a- a 
dx^ ' 
u 



t. 



stop 
V 



da d(j) 
E dx dx 



t. 



(11a) 
(lib) 
(He) 



stop 



The right hand side of equation (|llal) describes mass dif- 
fusion, £), induced by turbulent gas. The effects of this 
term have also been considered by Shariff & Cuzzi (per- 
sonal communication). The first term on the RHS of 
equations (jllbl) and (jllcp are drag forces on particle mo- 
tion relative to the static gas background. The aerody- 
namic stopping time, igtop is evaluated in terms of parti- 
cles sizes and gas densit ies in ^ 13.3.21 The last two terms 
on the RHS of equation (jllb[) are standard in analyses of 
gravitational instability (GI): the effective pressure ac- 
celeration from (isothermal) random particle velocities, 
c, and the acceleration from the disk's self-gravitational 
potential 4>. The validity of these model equations is 
discussed in §6.41 

We give the solution to Poisson's equation in terms of 
Fourier amplitudes which are denoted by tildes and have 
a spatiotemporal dependence oc exp(si -I- ikx): 



2TTGa 



T{kh) , 



(12) 



The softening term, T{kh) — 1/(1-1- fc/i) mimics the ef- 
fects of finite particle layer thickness, h ([Shu 1984). This 
factor connects the thin disk potential of long waves, 
kh <C 1, to the three dimensional solution for kh ^ 1. 
While the modes studied in this paper are generally long, 
T insures against the appearance of artificially short 
waves in a height-integrated model. 

3.1.1. Dimensionless Dispersion Relation 
We use the standard GI wavelength 

Ag = 2t:Ig = 2-K^GT,in^ (13) 

and the orbital time, 1//2, to define dimensionless veloc- 
ities, {f/, 1^} = {m, u}/(£g^), radial position X = x/Iq 
and time T = fit. The growth rate, s, and wavenumber, 
k = 27r/A, are non-dimensionalized as 

JC = k-£G^XG/>^, (14a) 
7 = s/J7, (14b) 

so the Fourier dependance becomes: exp(i/CA" +7T). In 
terms of dimensionless Fourier amplitudes, the equations 
of motion read 



+ K:'Qi,)a/S + ilCU = 0, 



(15a) 
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(7 + i/T,)i/~2y = I 

(7 + 1/t,)V + UI2 



1 — J^w o' 



(15b) 

0. (15c) 

The cubic dispersion relation for the (possibly complex) 
growth rate 7(/C) is 

(7 + r-'f + 1] (7 + IC'Ql) = (7 + r-i) (1 - J-w) . 

(16) 

We follow iWardI (|1976l 120001 ) in parameterizing the 
strength of self-gravity with the function 



^2TTGUkT{kh) + P 



(17a) 
(17b) 



Smaller J^w indicates stronger self gravity. 

At this point, the model requires four input parame- 
ters. These are the stopping time Ts plus three gravita- 
tional stability parameters; 



zf2 



ttGS 
ttGS 

ttGS 



0.2 



PR 



(18a) 
(18b) 

(18c) 



which incorporate the effective pressure (with T for 
Toomre), layer thickness (with R for Roche) and tur- 
bulent diffusion of solids (D), respectively. We express 
Qr in terms of the critical density of equation ^ and 

the midplane particle density p — E/(-\/27rft-). 

The three Q values are not independent. We relate 
them to the level of turbulent stirring in WS.2[ which has 
the added benefit of reducing the number of input pa- 
rameters to two. This approach allows the complete sur- 
vey of solutions in fJ5] First we describe some general 
properties of the dispersion relation. 

3.1.2. Standard GI: No Drag or Diffusion 

We recover standard axisymmetric GI by ignoring gas 
drag, Ts — >■ 00, and turbulent diffusion, Qu — 0. The 
dispersion relation takes the simple form 



(19) 

and a pair of 



7(7' + ^w) = , 

giving a neutral mode (NM) with 7 
density waves (DW) with 7 = ±V— -^w- In terms of the 
dimensional frequency lo = is, the DW solution obeys 
the familiar dispersion relation: 



2„2 



2TrGSk ■ T{kh) + k^c 



(20) 



The DW solutions correspond to traveling waves if J^w > 
0, and gravitational collapse for J- w < 0- In the th in disk 
limit where T — >■ 1 , we recover the lToomrd ^9M) insta- 
bility criterion, Qt < QT,crit = 1, with the marginally 
unstable wavenumber, /Ccrit = 1- In physical units 
A = Ag. 

We now consider the effects of our finite thickness cor- 
rection. First note that T enforces a Roche-like limit. 
Even for a perfectly cold layer with Qt = 0, instability 



{J-w < 0) requires Qr < 2. If we fix Qt = Qr (the 
standard c = Qh relation from Keplerian orbital oscil- 
lations or hydrostatic balance) then Qt.ciH ~ 0.55 and 
/C(,j.ij sa 1.2. This solution ignores the compression of the 
layer by self-gravity which enhances the vertical gravity 
by a factor 



27rG'I] 2 



(21) 



Vertical hydrostatic balance then gives a relation Qt = 
<3r\/1 + 2/Qr (as in ImTdl [200(1 . In this case 
Qx.crit ~ 0.77 with QR.crit ~ 0.26 and /Ccrit ~ 
1.05, in good agreement with more rigorous analyses 
(|Goldreich fc Lvnden-B"elilll965D . Thus the finite thick- 
ness correction to Toomre instability is relatively mod- 
est when vertical compression is included. Secular GI 
has much longer wavelengths, making the finite thick- 
ness corrections insignificant in most cases of interest. 

3.1.3. Proof of Generic Instability 

We now prove that the full dispersion relation of equa- 
tion (jl6p always yields a growing mode for some wave- 
lengths, for any choice of stability parameters or stop- 
ping time. This result is fundamentally different from 
standard GI, as we just reviewed. The existence of a 
growing mode does not assure its astrophysical relevance, 
so we impose extra constraints in ij4.1l Nevertheless, the 
following proof demonstrates that instability is quite ro- 
bust and does not de pen d, e.g. on the specific turbulence 
model we adopt in §3.21 Moreover we know that a nu- 
merical root finding algorithm has failed if it does not 
identify a growing mode. 

Since the dispersion relation is a cubic in 7 with real 
coefficients, we can use Descartes' "rule of signs" to show 
that there is always at least one real, positive root (and 
thus a growing mode). We write equation (|16|) as 



7^ + C27' + Ci7 + Co = (22) 



with 



C2 = 2/ts + /C2q2^ > 

Co = (-Fw-l)/r, + C'/CVr, 

We use the long wave limit where J^w 1— 2/C as /C ^ 0. 
We have C2 > and Ci > since > for long 
waves. The Co coefficient is negative if C'K? < 1 — J^W: 
which holds because G'K? < 2JC as /C — > 0. Thus we can 
ensure (for some JC) one, and only one, sign change in the 
coefficients, which guarantees a positive root. Q.E.D. 

3.2. Turbulence Model 

We will now describe the response of solid particles 
to turbulent stirring. This section culminates in equa- 
tion p3l) , which expresses the self-gravity parameters 
Qt, Qr and Qd in terms of the stopping time Tg and 
a parameter 



rf2 



nGS 



ttGS 



(23) 
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that compares turbulent stirring to particle self-gravity. 
Here Dg is the strength of turbulent diffusion in the gas, 
and 

a = Dgf2/cl (24) 

is a more familiar dimensionless measure of turbulent 
intensity with Cg the sound speed in the gas. We discuss 
in 33 why this diffusive a can be much smaller than the 
(purposefully similar, but not equivalent) parameter used 
in studies of accretion disks. 

This subsection makes heavy use of the results of 
lYoudin fc Lithwickl (|2007l hereafter YL07). We extend 
their results in our equations (f29|) -- ([3T|) by including self- 
gravity in the scale height calculation. 

3.2.1. Turbulence and a Values 

We adopt a relatively simple model for disk turbu- 
lence. As in Kolmogorov turbulence, the largest eddies 
dominate the energy budget. We further assume that 
icddy = l/-!^, i-G. the turnover time of these integral scale 
eddies equals the orbital time scale. This assumption re- 
quires a turbulent speed T4ddy — V^'^g ^-^^d eddy length 
4ddy = y^Cg/n SO that toddy = ^ddy/Voddy and the 
gas diffusivity = T 4ddv^cddv satisfies equation (j24p . as 
noted by Cuzz i et all (jloOl). 

We further assume that the turbulent kinetic energy 
is isotropic with partial correlation between the radial, 
Ug, and azimuthal, Ug, turbulent speeds. Specifically, we 
take 

(Ug) = (vg) = (wg) = 4(ugt;g) = Foddy = aCg (25) 

where Wg is the vertical turbulent speed and angular 
brackets indicate time-averages. While the velocity cor- 
relations are of modest significance for particle stirring 
(see YL07) they are a primary concern for the angular 
momentum transport and the accretion of gas onto the 
protostar, see SJT] 

3.2.2. Particle Response to Turbulence 

The radial particle diffusion coefficient is given by 
equation (36) of YL07 . With our equation (1^ 



D 



1 



-4r2 



(1 



(26) 



The detailed Ts dependence will vary with the exact prop- 
erties of the turbulent cascade. The important result is 
that Tg <C 1 solids diffuse at a rate similar to any passive 
contaminant, D ~ Dg. Particle clumping in turbulence 
could lower D, but we do not have a good model for 
this effect. As particles decouple and Tg ^ 1, diffusivity 
drops as D/Dg cx t~'^. Thus a rough value of the so 
called Schmidt number is Sc = Dg/D ^ 1 + t^/A. This 
result differs from the often us e d, bu t incorrect, value of 
SccDC = 1 + Ts in lCuzzi et all (|1993D . see YL07. 

Since we are interested in radial pressure support, the 
relevant particle dispersion speed, c, is the radial com- 
ponent. We use 



+ 2t2 + (5/4)r3 

c = — — ^ VaCff 

1 + T/ 



(27) 



A similar result for c follows from equation (33a) of YL07 
and our equation (j25p . In equation (j27p we adopt a sim- 
pler Ts dependence. This choice is purely for convenience 



and gives indistingui shable num e rical r esults. These de- 
tails are not crucial. iVolk et all ()1980D showed that the 
same basic scaling, c ~ Voddy/(l + 7's)^^^, holds when the 
orbital dynamics included by YL07 are ignored. 

The particle scale height, h, usually obeys the well 
known result 



h — hstd = 



(28) 



where hg = Cg/ f2 is the hydrostatic gas scale height (for 
weak self-gravity). Intuitively the particle layer narrows 
for weaker t urbul ence or more loosely coupled solids. 
iCuzzi et"aLl ()1993f ) derive the result for Tg ^ 1 by com- 
paring settling and diffusive time scales. iCarballido et al.l 
(|2006) show that the same scaling also holds for Tg ^ 1 
using both simulations and velocity diffusion arguments. 

Corrections to hstd arise from self-gravity and the per- 
fect mixing limit. Including self- gravity gives 



h = — 




(29) 



with the amplification of vertical gravity by ip from equa- 
tion (|2ip . To derive this result (which as far as we know 
is new) for Tg <C 1, equate the shortened settling time 
isctt — l/(^TgV'') with the diffusion time tdis = h? /Dg. 
The possibility that turbulent intensity varies with self- 
gravity is not considered here, but could be included to 
model gravito-turbulence. The same expression holds for 
Tg ^ 1, where the vertical oscillation frequency shortens 
to ujz — fi/\fip. The random velocity, Cz — ya/r^Cg, 
imparted by turbulent fluctuations gives h^p = Cz/uiz, 
which obeys equation ([^ . Since i/' itself depends on h, 
the explicit solution reads 



a 

Ts 



e% + -hi 



(30) 



which reduces to equation ([28)) for small £g : which is de- 
fined in equation (jl3p . 

Perfect mixing with the gas imposes an upper limit 
h < hg. Following YL07 [their equation (28)] we enforce 
this limit to give our most general result 



^full 



TgV" + a 



(31) 



Again there is an implicit h dependence in ip. 

The ratio of particle to gas density, important for drag 
feedback effects, is 



p„ Tia h \ a 



(32) 



which uses hstd for illustration, though equation (1511) is 
used in numerical calculations. 

3.2.3. Dimensionless Stability Parameters 

We can now express the Q parameters from equa- 
tion (jl8p in terms of Qa and Tg. Applying equations (|26p. 
((27)) and ([31]) gives 



v/1 + 2rg2 + (5/4)7 



(33a) 
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v/l + Q2/r,(l + a/T,)-l 



1 + a/i 



4r2 



1 + t2 



(33b) 



(33c) 



Note that equation (I33b[) depends on a in addition to 
Qa and Ts. To reduce to only two parameters, we take 
the a ^ Tg hmit, equivalent to using equation (|29l) . and 
obtain 

Qr « Vl + Qlh. - 1 . (34) 

This approximation is acceptable on a few counts. First 
it is conservative, since we are removing the restriction 
that h < hg and allowing artificially low particle densi- 
ties. Second, relevant growth rates typically have a < Ts, 
as seen in Figs. [Hand HI Finally, even for a > Tg, wave- 
lengths are long enough that the finite thickness correc- 
tion is negligible, kh — /CQr ^ 1, even with the artifi- 
cially enhanced h. Thus equation (|34p is a lway s accept- 
able for our study, even though equation (j33bl) is more 
general. 

For Tg <C f , Qt — Qd — Qa is particularly simple. It 
reflects the result that tightly coupled particles share the 
same velocity and diffusion coefficient as the dominant 
turbulent eddies. Furtherore Qr — Qal \f% ^ 1 for 
weak self-gravity, reflecting the fact that h ^ cQ for 
turbulently-stirred small particles. 

3.3. Disk Model 
3.3.1. A Scalable Minimum Mass Nebula 

For numerical evaluations we adopt a reference mini- 
mum mass Solar Nebula model, with scaling factors to 
explore the effect of changing disk parameters. In our for- 
mulation (similar to CYIO which has further discussion) 
the surface density of solids and gas and the midplane 
temperature are 



T 



20 Z%F Rll^' g/cm^ 
2000Fi?lf/^ g/cm^, 



"•AU 



(35a) 
(35b) 
(35c) 



where i?AU = i?/AU is the disk radius in AU. The 
"metallicity" Z = S/Sg = 0.01Z% and the mass nor- 
malization F are two important adjustable parameters. 
Water and methane condensation at "ice lines" can give 
sharp transitions in E. We do not model this in detail, 
because a range of enrichment processes also affect Z. 

We adopt a cool disk temperature characteristic of pas- 
sive flared disks around Sun-like stars, and fix = 1 
in this work. Our profile is modestly steeper than the 
T cx of (jChiang fc GoldreichI [T997[) . just to give 

analytic s calings that are both simpler and similar to 
the classic iHavashil (|198lD model. This idealization has 
little effect, given the weak dependence of our results on 
such small temperature changes. 

The disk mass within a cutoff radius -Rout is 

Mdisk = 0.02 F(l + Z) v/i?out/(50 AU) Mq. (36) 

The gas disk is gravitationally stable since its Toomre 
parameter. 



(37) 



satisfies Qg > 1, with rn* the stellar mass in solar masses. 
Our measure of particle self-gravity relative to turbulent 
stirring, from equation (j23p . 



1.3 



10- 



FZ: 



% 



AU 



exceeds unity in most disk models. 
The midplane gas density is 



Pg 



27r/i„ 



V /t cm3 



(38) 



(39) 



for a vertically isothermal disk, ignoring compres- 
sion from the weight or self-gravity of t he solids 
(|Nakagawa et al.lfT986t IWakita fc Sekivall2008D . 
The mean free path for molecular collisions, 



W^F^AU cm, 



(40) 



is inversely proportional to pg and important for the drag 
forces that we now consider. 

3.3.2. Particle Properties and Gas Drag 

Our model considers a population of particles with 
uniform radius, a, and internal density, — 2 g/cm"^. 
The general dependence of drag forces on particle, gas 
and flow properties is complex, but appro ximate formu- 
lae apply to distinct parameter regimes (jAdachi et al.l 
1976; WcidcnschiUing 1977). The stopping time igtop = 
AV/ /drag is the ratio of the relative flow speed to the drag 
acceleration. Fig. [1] plots the dimensionless Tg = ^?tstop 
for particle sizes from a = 1 mm to 30 cm for both a 
minimum mass disk with F = 1 and a gas depleted disk 
with F = 0.1. 

Epstein drag is the most relevant regime. It applies for 
small sizes, a < (9/4)i?g, and low gas densities, with 



TEp = 



p,a 



2ttX;^ 



2.5 X 10 



-40111111 p3/2 



F 



R"^, (41) 



where = a/mm. Epstein drag is ready identifiable 
in Fig. [T] as the region where Ts increases with R (outside 
0.3 AU for 1 mm and outside ~ 2 AU for 30 cm). 

Stokes drag applies to particles too large for Epstein 
drag, and a Reynolds number, Re = a/SV / {cg^g) < 1. 
The stopping time for viscous Stokes drag 



TSto = TEp 



4a 
94 



1 X 10 




(42) 



applies in the inner disk, where gas densities are high, 
and for small solids that do not experience a large head- 
wind. In Fig. [T] we identify the Stokes regime for 1 mm 
and 3 cm solids where Tg declines sharply with radius. 
The brief overlap of the F = 1.0 and F = 0.1 (gas de- 
pleted) curves demonstrates that rsto does not depend 
on gas density. 

In Fig.lll turbulent drag applies in regions where Tg has 
a relatively flat dependence on R. Appendix [A1 discusses 
non-linear drag regimes. Unless stated otherwise, the 
numerical estimates and scalings in this work use Epstein 
drag. 
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3.3.3. Radtal Drift 

Particles experience a net inward drift because they 
encounter a headwind in pressure-supported gas disks. 
The orbital speed of a purely gas disk is slower than 
Keplerian by a speed rjQR, where 



■n ■ 



dP/dR 
2pgf2^R 



1.3 X lO-^izl/Rl^ (43) 



measures the radial gradient of the gas pressure, P. 

The radial drift time scale, defined as the time to reach 
the star if the speed were to remain constant, is 



^drift = 



R 



\Vdnft\ 



= /i: 



1 



(44) 



This equation includes the inertial correctio n factor 
(|Nakagawa et al.lflQSl lYoudin fc Chianell2004[ ) 



/i. 



(ptot/Pef = (1 + p/pgf 



(45) 



which is unity in the test particle limit that, p <C Pg 
Test particles with Tg = 1 have the shortest drift time 



min(t 



drift ) 



120v^i?Au//T yr. (46) 



The "meter-size barrier" refers roughly to the size of Tg = 
1 bodies in the inner disk. As Fig. [T] shows, ~ cm-sized 
bodies drift fastest in the outer disk. Particle inertia 
both lengthens the shortest drift time — as min(tdrift) oc 
VTin = Ptot/pg — and shifts it to larger Tg = ptot/pg. 
In the small particle particle regime, with Tg < 1 



idr, tight 



fn 



2rjQTs 



2.5 X 10'^ 



•AU 



(47) 



Larger particles drift faster in the test particle limit. 
However inertia can cancel this trend. With a fixed level 
of weak turbulence that ensures p ^ pg, equation p2p 
gives /in « Z'^Ts/a, canceling the explicit size depen- 
dence above. 

For larger particles with Tg > 1 , the drift time increases 
with size and radial distance as 



^dr. loose 



2.7 X 10^ 



(48) 



where R^o = i?/(50AU) and = a/(lm). There is 
no inertial correction for large solids with Tg Ptot/Pg, 
because drag forces are too weak to slow the headwind. 

4. NUMERICAL RESULTS 

We now present the main numerical results for the 
gravitational collapse of solids in a turbulent gas disk. 
The fiducial protoplanetary disk models of !j3.3l provide 
the physical values. Our most commonly plotted quanti- 
ties are the growth time, wavelength and mass of solids 
in a collapsing ring, which we define as 



M ■ 



grow 

A 



[i7Real(7)]-i , 
2ti\RY. . 



(49a) 
(49b) 
(49c) 



To obtain these quantities, we solve the dispersion re- 
lation in equation (jl6p for the fastest growing modes. 
This solution involves maximizing the real part of the 



a = 1 mm 



10*^ 



_ 10'' 
I 10' 



10^ 
1 

10^ 
1 

10-^ 














10 p 


tr 


1 r 


to 


LU 


10-% 




10-^ r 



lo-'^ 

10 



R/{2H,) 





1 10 

R[AU] 



Figure 2. Properties of the dissipative gravitational collapse of 
mm-sized particles vs. disk radius in the reference model for several 
values of the turbulent diffusion parameter a = 10"^, 10~*, 10"^'' 
{solid, dashed, and dot-dashed curves). Grey sections of curves 
denotes when growth is too slow. While remarkably low, these 
a values do not preclude stronger turbulence away from the mid- 
plane, and also ignore the possibility that small-scale clumping 
reduces diffusion. Top: Growth times in years. Middle: Wave- 
lengths in AU compared to (half) the disk radius, the gas scale 
height kg, and the standard GI wavelength \q (brown, black and 
grey dotted lines, respectively). Bottom panel: Mass of collaps- 
ing ring in Earth masses. Subsequent fragm entation yeilds much 
smaller planetesimal masses. See text ( >]4.2II for details. 



(dimensionless) growth rate 7 with respect to the (dimen- 
sionless) wavenumber JC, see Sj5]for details. We quote a 
ring mass because the modes are axisymmetric, but em- 
phasize that this is not a predicted planetesimal mass. 
Gravitational fragmentation during later stages of col- 
lapse would likely produce smaller bodies. 
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a= 10"^ 




1 10 
R[AU] 



10^ 



Figure 3. Similar to Fig.[2]with turbulence held fixed at a = 10~® 
and other parameters varied. The reference model with a = 1 mm 
solids [solid black curve) is shown again for comparison. Larger 
a = 3 cm solids {green dot-dashed curves) give faster and smaller 
scale growth. Increasing the particle-to-gas surface density ratio to 
Z = 0.1 also promotes growth, whether the gas content stays fixed 
{F = 1.0, red dotted curve) or is depleted (F = 0.1, blue dashed 
curve). With both gas depletion and 3 cm particles {purple triple- 
dot-dashed curve) growth is even faster and smaller scale. The 
kink in A and M^i^g at 20 AU is a real feature (not a numerical 
artifact) as explained in i)5.4l 



4.1. When is Collapse Relevant? 

Since our instability always produces growth (see 
t j3.1.3|) . we require the following conditions be met for 
collapse to be relevant: 



idrift : 
ilifc = 

R/2. 



10^ yrs , 



(50a) 
(50b) 
(50c) 



Growth must be faster than both the radial drift time 
of equation (|44)) and an assumed 10^ year disk lifetime. 
Also the wavelength must fit comfortably inside the disk. 
The drift requirement is usually the most stringent. 

The characteristic growth time, tgrow, only represents 
a single e-folding. Many e-foldings are generally required 
for significant density growth. The actual time to form 
planetesimals depends on seed perturbation amplitudes 
and non-linear evolution, both beyond this analysis. To 
avoid introducing another uncertain parameter, we sim- 
ply comment that many e-foldings are possible because 
growth could proceed over several Adrift, especially since 
the growth of particle overdensities will slow the actual 
drift timescales. Moreover disk lifetimes could extend 
to ~ 10 Myr, though this is less often the relevant con- 
straint. These issues touch on the application of a local 
instability to a globally evolving disk, as discussed fur- 
ther in gOl 

4.2. Growth for Fixed a 

Our reference model consists of particles with radius 
a = 1 mm in our minimum mass disk model with _F = 1 
and Z = 0.01. Fig.[2]plots growth times and wavelengths 
for different levels of the turbulent parameter a with all 
other disk quantities fixed. Curves become light grey 
whenever the conditions of equation (|50p are not met. 
Low a values are required for collapse in our nominal 
disk. With a = 10"^ collapse only proceeds outside 10 
AU. For a = 10~* collapse occurs everywhere except a 
region near R — 0.3 AU, where the peak in growth times 
corresponds with the minima in stopping times at the 
Epstein-Stokes transition, see Fig. [TJ 

Increasing the strength of turbulence gives slower col- 
lapse that acts over long wavelengths. Our physical 
derivation of equations ^ and (fTO|) explains the basic 
scalings. Including the more precise coefficients of equa- 
tion (|56|) , the growth time for Tg ^ 1 becomes 



1 aC 



(51) 



The more rapid collapse in the outer disk occurs mainly 
due to weaker drag coupling, with larger Tg. This aero- 
dynamic effect overcomes the radial increase in orbital 
time scales. Stronger self-gravity, as measured by Qg, 
also assists growth in the outer disk. 

Substituting numerical scalings for Epstein drag and 
the reference disk, the growth properties are 



A 

AU 

Mring 



-R 



4 X lO"" 



Zo/n„ 



AU ■ 



Mff 



0.13 



%"mm 
Omm-RAU 



(52a) 
(52b) 
(52c) 



where turbulence is normalized to a_8 = a/10~^. Dif- 
ferent scalings (not given) apply to the Stokes regime at 
small R. 

Significantly faster and smaller scale growth is possi- 
ble with deviations from the reference model. Fig. [3] ex- 
plores the effect of varying disk parameters while keeping 
a = 10~® fixed. Larger particle radii with bigger Ts give 
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a = a„ 



1 

10-^ 



fl = 30 cm & gas dep. 
a = 3 cm & gas dep. 
a = 30 cm 
a = 3 cm 




10" 



10 



10^ 



R[A\J] 



Figure 4. The maximum strength of turbulence, given by Qmaxi 
which allows collapse. T he c riteria for sufficiently rapid collapse 
are described in equation II50I I. The Omax values increase with par- 
ticle size as shown for a = 1 mm {black solid curve), a = 3 cm 
{dash-dotted green curve) and a = 30 cm {dotted red curve). The 
enrichment of the solid to gas ratio also increases Qniax- We show 
this for gas depleted models with Z = 0.1 and F = O.l for both 
a = 3 cm {triple-dot dashed purple curve) and a = 30 cm {dashed 
blue curve). See Fig.[5]for the properties of modes with a = Omax- 



much faster growth, as the a — 3 cm (^dot-dashed green) 
curves demonstrate. Deviations from powerlaw behavior 
at large R arise from the transition to loose coupling, 
Ts > 1, which f}5] explores in more detail. 

Higher particle fractions, Z, produce faster growth 
by increasing particle self-gravity. As an added bene- 
fit, more particle inertia slows radial drift. Whether Z 
increases by adding solids or depleting gas (or some com- 
bination) the effect on growth rates and w avelengths is 
similar. Indeed equations (|52ap and (|52bp are indepen- 
dent of F. While a higher mass disk has stronger particle 
self-gravity (at fixed Z), higher gas densities cancel this 
effect by giving tighter coupling (in the Epstein, but not 
Stokes, regime). The ring mass in equation (I52cp does 
increase with F. At fixed Z, there is obviously more 
particle mass in a fixed area when the total mass rises. 

Fig. [3] demonstrates these metallicity effects for two 
high Z ~ 0.1 models, with F — 1 (red dotted curve) and 
F = 0.1 (blue dashed curve). The overlap in tgrow and 
A is evident at intermediate R. The curves separate at 
small R with the transition to Stokes drag, and at large 
R with the transition to Tg > 1, which occurs at smaller 
R in the gas depleted model. 

The effect of particle growth and enhanced metallic- 
ity can interact constructively. The purple (triple dot- 
dashed) curve labelled "all" shows this result by com- 
bining the a = 3 cm, Z = 0.1 and F = 0.1 parameters. 
Growth times are even faster and ring masses are much 
smaller when these effects are considered together. 

4.3. Maximum Allowed Turbulence 




a = 30 cm & gas dep. 

= 3 cm & gas dep. 
fl = 30 cm 
fl = 3 cm 
fl = 1 mm 



10^ 



R[AU] 



Figure 5. Properties of modes with the maximum allowed turbu- 
lence, a = Cfmax, for the same models as in Fig. [4] We can visually 
identify which constraint sets amax '■ if t grow — 10^ yr, it is the 
disk lifetime, and when A = R/2, it is the disk size. Other wise the 
drift time scale is the relevant constraint. See the text ( i|4.3l l for 
an explanation of the overlapping wavelengths and masses. 



Instead of fixing the level of turbulence, we now de- 
termine the strongest level of turbulence that allows col- 
lapse. Fig. 2] plots aniax; the largest a value for which 
growing modes satisfy all the constraints of equation ([SO)) . 
Not surprisingly, the same factors that give faster growth 
— larger solids and higher particle mass fractions — also 
yield larger amax- 

The constraints on turbulence are strongest inside a 
few AU, where a < 10~^ is required for sizes up to 30 cm, 
even with some metallicity enrichment. While the outer 
disk generally allows collapse with stronger turbulence, 
it still requires a < 10"'^. 

Analytic fits to amax come with even more caveats 
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a = a„ 



DC 

Q. 



-5- 10h - - 

Q. 




1 

10-^ 



1 10 10^ 

i?[AU] 



Figure 6. The minimum particle density, p, that gives collapse, 
obtained by setting a = Qmax, using the same disk models as 
in Figs. HI and [51 (Top): p plotted relative to the Roche density 
[equation Q], demonstrating that secular GI occurs well below this 
standard threshold. (Bottom): p plotted relative to gas density to 
show that drag feedback is always significant. 

than usual, but are useful for a rough understanding. 
When drift speeds are the limiting constraint, the condi- 
tion t„row < idrift SCtS 



^sZ^ fin 

2vQl 



0.6 X 10 



-8' 



Pt 



R 



3/2 
AU ■ 



(53) 



for Tg <C 1. (Again the numerical scalings apply only to 
Epstein drag.) This basic scaling explains many of the 
curves in Fig. 21 e.g. the a = 1 mm curve beyond 1 AU 
and segments of the a = 3 cm curves. For transitional 
Ts ~ 1 values, amax has a flat R dependence, e.g. for the 
a = 3 cm curve beyond 10 AU. For Ts ^ 1, a similar con- 
straint to equation (15^ holds, again with ctmax (x R^/"^ 
as seen in the a = 30 cm curve beyond 20 AU. The kink 
near 60 AU in the gas depleted a = 30 cm curve occurs 
because the disk lifetime constraint becomes relevant. 

The upper limit to allowed turbulence, amax, varies 
greatly with disk model and has no absolute upper limit. 
The range of models we consider all require amax ^ 10"'^, 
and lower values in the inner disk. 

4.3.1. Mode Properties 

To better understand these limits on turbulence, we 
plot the properties of modes with a — amax in Fig. [5l 



These growth times, wavelengths and masses here take 
the largest allowed values for a given disk model. For 
a < amax, collapse is more rapid and smaller scale. 

Fig. [5] allows us to readily identify the relevant con- 
straint that sets amax- The 10^ year lifetime constraint 
applies when the drift time scale is even longer. As shown 
in the top panel, this occurs both for loose coupling (be- 
yond 60 AU in the gas depleted a = 30 cm model) and 
for very tight coupling (mm-sized solids near ~ 0.3 AU 
at the Stokes-Epstein transition). In the middle panel, 
we see that the disk size constrains wavelengths in our 
enhanced metallicity models, especially in the inner disk. 
The radial drift time scale is the relevant constraint on 
amax wherever igrow 7^ 10^ yr and A ^ i?/2. 

A striking feature of Fig.[5]is the clustering of the wave- 
length (and also mass) curves for modes for models with 
very different Tg. Substitution of amax oc Tg from equa- 
tion (|53p into the wavelength result of equation ()56|) gives 

A(amax) « « 0.017 ^'"^^^°^^^ AU , (54) 

which is independent of Tg. Thus the maximum wave- 
length imposed by drift is at least l/(2?7) larger than the 
standard GI wavelength, and more if the inertial correc- 
tion factor, /in, is significant. The gas depleted curves in 
Fig. [5] have larger maximum wavelengths when Tg < 1. 
In these cases, higher particle inertia increases /in values, 
raising A(aniax) to the A = R/2 limit. 

4.3.2. Minimum Particle Density 

The volume density of solids, p, is both a measure of 
self-gravity — though an indirect and imprecise one for 
secular GI — and a measure of feedback effects that give 
rise to drag instabilities. When a — amax we find the 
minimum p required for collapse. 

The top panel of Fig. [B] plots p/pR — 0.2/(3r. Since 
P ^ PR, secular GI occurs well below this standard 
threshold. While !j3.1.3l proved that collapse formally 
occurs for any and thus p, these numerical results 
demonstrate collapse for small p when astrophysical con- 
straints are imposed. Below the p values plotted in Fig.|6l 
collapse is too slow to be relevant. 

We emphasize that density is an indirect measure of 
the relevance of secular GI. Radial diffusion is the dom- 
inant stabilizing influence, so the vertical diffusion that 
sets p is more of a side effect. If we consider asymmetric 
turbulence and vary just the vertical diffusion, p would 
change, but growth rates would be unaffected because 
the finite thickness correction is negligible for long wave- 
length modes. Drift rates however would change from 
the inertial correction. 

The bottom panel of Fig. [5] plots the particle density 
— again, the minimum value required for collapse — 
relative to the gas density. Remarkably, the minimum 
mass models with = 1.0 all have p ^ Pg- Applying the 
drift constraint of equation (|53p to the density relation in 
equation ([5^ gives 



p(an 



2.0- 



(55) 



Thus p/pg is unity for _F 1 models (when Tg < 1) and 
larger for the gas depleted F = 0.1 models, explaining 
the bottom panel Fig. |6l It is a numerical coincidence 
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Table 1 

Growth across parameter space 



# 


Drag 


Turbulence 


I 


Ts < 1 




II 


Ts > 1 


Qa > 


III 


Ts > 1 




IV 


Ts > 1 




V 


Ts < 1 


^ ^ 3/2 



Tf/Q^ 

0.2/Ql 



0.6 



Ts/Q^ 

0.2ts/Q2 



0.4 



-2/3 



4/3 



Note. — Col. (1): Region label, see also Fig.H Cols. (2) 
& (3): Region boundaries in Tb-Qc space. Cols. (4) & (5): 
Growth rate and wavenumber of fastest growing modes. 

— unrelated to the mechanism of the instability — that 
standard minimum mass gas disks give Qg^/Sry of order 
unity. But this coincidence is significant. The conse- 
quence that p ^ Pg for gravitational collapse means that 
drag feedback is relevant, and can drive particle-driven 
turbulence. See ^6.31 for more discussion. 



5. GROWTH ACROSS PARAMETER SPACE 

We now survey the behavior of our instabilities across 
the full range of parameter space. Varying only two pa- 
rameters — the stopping time, Tg, and the strength of 
turbulence relative to particle self-gravity, Qa — covers 
all possibilities, thanks to the turbulence model of ^3.2[ 
Once solutions to the dispersion relation of equation (|16p 
are tabulated, translating the conditions in an arbitrary 
disk model to appropriate values of Tg and Qa allows sim- 
ple look-up of the solutions. Furthermore, we gain both 
a deeper understanding and a useful check on numeri- 
cal solutions by understanding the analytic behavior for 
extreme values of the input parameters. 

Fig. [7] maps the growth rate, 7fgm, and wavenumber, 
A^fgm, of the fastest growing modes versus Qa and Tg. 
We obtain these solutions by maximizing each of the 
three roots of equation (IT5|) with respect to /C over a well- 
sampled grid of Ts and Qa valuesEI 

Fig. [7] shows the general trend that both growth rates 
and wavenumbers tend to increase for small Qa (weaker 
turbulence or stronger self-gravity) and larger Tg (bigger 
solids or weaker drag). However the behavior has a differ- 
ent character in different regions of parameter space. We 
identify five distinct regions of parameter space, shown 
in Fig. m Table [1] lists the analytic growth rates and 
wavenumbers by region. 

We now describe the behavior by regions, with less 
relevant cases deferred to appendix |B] The transition 
from secular to standard GI occurs across the region H 
to HI boundary. Fig. [7] shows the corresponding sharp 
rise in growth rates, explored further in ^15. 41 

5.1. Region I: Secular, Tight 
Region I is the most relevant for planetesimal forma- 

3 /2 

tion. With Tg ^ 1 and Qq » Tg , it describes tight cou- 
pling and turbulence that is not incredibly weak|3 The 

3 For extreme values of the input parameters, renormalizing K 
towards a known solution as K. ^ K.' = /CQ^/ts helps avoid nu- 
merical difficulties. 

* Appendix IB. II shows that the weak turbulence solutions in 
region V are not relevant. 



physical derivation of secular GI in ^ — a simple balance 
between terminal velocity infall and diffusion — applies 
in region I. 

The fastest growth rate and corresponding wavenum- 
ber in region I, 



(56) 



are limiting solutions of the full dispersion relation. They 
are more exact versions of the physically derived equa- 
tions (IH) and (fTO| . The smooth variation with Qa, also 
evident in Fig. [71 demonstrates the insignificance of the 
Qa = Qt — 1 transition for Tg ^ 1. 

Growth rates in region I can be as fast as 7 ~ 1 /rg ^ 1 
at the region V transition. Nevertheless, growth rates are 
slower than dynamical, 7^1, because our disk models 
have Qa > 1 as in equation (|38l). 

The thin disk approximation holds in region I. From 
equations and (1551) . /CiQr, ^ 1 for weak self-gravity 
with Qa > ^/t\- The strong self-gravity case of Qa < 
^ 1, is less relevant, but still gives only a modest 
finite thickness correction, /CiQr ~ 1. (And thankfully 
we need not sub-divide region I.) 

Density waves are rapidly damped by drag in region 
I. Growth occurs via the neutral mode, introduced in 
ij3.1.2l as trivial in the absence of drag. The eigenfunc- 
tions of the particle response are Ui = 2Vi/ts = 2zTgO'i. 
This strong density response, with a ^ U ^ V, makes 
it unlikely that secondary flow instabilities will halt the 
gravitational collapse once it begins. 

Wavenumbers other than the fastest growing can also 
collapse, at a rate determined by the leading order dis- 
persion relation: 



71 (/C) « 2rg/C-/C2 



(57) 



This result has implications for how collapse proceeds for 
a distribution of particle sizes as discussed in 



5.2. Region II: Secular, Loose 

Region II describes larger, loosely coupled sohds with 
Tg 3> 1 and relatively weak self-gravity characterized by 
Toomre's Qt — Qa/^/r^ ^ 1- Like the tight coupling 
case of region I, gravitational collapse involves the secular 
destabilization of a neutral mode. 

The physical derivation of region II growth follows 
the qua htative description o f secula r GI in the introduc- 
tion of [Goodm^XPindo3 dlOOO) ■ As in fj2l we con- 
sider a ringlike perturbation of the solids with radial 
self-gravity gr ~ Ga. The orbital speed at the outer (in- 
ner) edge of the ring increases (decreases, respectively) 
by w ~ gr/{2fi) in response. Assuming the gas is unaf- 
fected by this acceleration (see §6.5|) . the extra headwind 
(tailwind for the inner edge) causes a radial drift towards 
the ring center. By angular momentum conservation this 
drift speed is 



2v 



Ga 



stop 



(58) 



Due to the weaker coupling this scaling differs from equa- 
tion (jlj . The collapse rate again follows from mass con- 
tinuity as 
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Figure 7. {Solid curves and grayscale): Contours of (the log of) the dimensionless growth rate, 7, {left panel) and wavenumber, JC, {right 
panel) of the fastest growing mode as a function of the strength of the turbulence - measured by Qa - and the size of particles - measured 
by Ts. The dotted curve denotes Qt = Ij the nominal threshold for standard GI. Slow growth occurs for Qt > 1 (above the dotted curve) 
and the Qt = 1 contour only has significance for rs 2> 1 • See Table [T] and ij5]for explanations of different parameter regimes. 



As in f}2l we consider the minimum A allowed by pres- 
sure and diffusive stabilization. It turns out that both 
are of roughly equal importance for loose coupling (and 
our turbulent stirring model). Pressure stabilization sets 



Larger particles 



a minimum wavelength 




1 

= at. 



stop 



Figure 8. Overview of Qa vs. Ts parameter space. Region labels I 
— V are used in Tableland the text. Regions can be identified in 
Fig.[7l which is not centered on (1,1). The solid line that divides 
region I and runs between regions II and III) denotes Q-p = 1, 
marginal stability for standard GI. 



A, 



c 



"eddy 



(60) 



where we use the result c ^ V^ddy/ from equation (|27p 
for the random speed of Tg ^ 1 solids stirred by turbu- 
lence. Mass diffusion imposes a time scale constraint 
sX^ /Dp > 1 that also sets a minimum wavelength 



Xd 



f2Dg 
GSt, 



GSt, 



(61) 



using D ^ Dg/r^ for the turbulent diffusion of large 
bodies, see equation (f26l) . The agreement with equa- 
tion (|60|) shows that pressure and diffusive stabilization 
are of roughly equal importance. 

The resulting growth rate,sii ^ fllQ\, is independent 
of stopping time. As Tg increases, shorter wavelengths 
collapse, but the rate is unchanged because there is dis- 
sipation to drive infall. In dimensionless units the fastest 
growing modes in region II obey 



711 = ^ii/t-s = Cfc/Q 



2 

QL ' 



(62) 



where the numerical coefhcient Cfe ~ 4/21 reflects the 
combined effects of pressure and diffusive stabilization 
from an expansion of the full dispersion relation. Aside 
from this numerical coefhcient, the wavenumber has the 
same basic scahng K ^ t^IQ\ in regions I and II. 

In region II, growth rates are always slower than dy- 
namical with 7n <C 1. Furthermore, 7iirs ^ 1, so 
that drag can dissipate excess energy generated dur- 
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ing collapse. The thin disk approximation holds with 
ICiiQr - y^/Qa - 1/Qt < 1. 

5.3. Region III: Recovery of Standard GI 

Region III is where our model reproduces standard GI 
due to unstable density waves as in ij3.1.2l Radial mass 
diffusion is negligible and drag forces (still relevant in 
region II despite Ts ^ 1) are weak compared to particle 
inertia. Turbulence still sets velocity dispersions (giving 
Qt) and the particle layer thickness (giving Qr). The 
boundaries of region III are 1 / y/r^ <C Qa ^ \/^' "which 
requires Ts ^ 1. 

Since Qt ~ Qa/y/r^ ^ 1, region III is likely not acces- 
sible in practi ce. Self-gravitating disks l i kely self-regulate 
to Qt J 1 (IGoldreich fc Lvnden-Belll [19651: ISkilllQQl 
lGanimiell2001[ ). Thus the boundarv between regions II 
and III (described below) is likely more relevant than the 
heart of region III. Unlike reg ions I and II, G I might not 
result in collapse in region IIL IGammid (|2001[) shows that 
energy dissipation must occur on orbital time scales for 
standard GI to give collapse. With Tg ^ 1, dissipation 
by gas drag is not sufhciently rapid. 

Linear growth in region III occurs at the free-fall rate 



(63) 



The minimum wavelength is set by effective pressure sup- 
port, A < Ap, from equation (j60l) . which gives a growth 
rate sm ~ Q/Q^. Diffusive support is too small scale. 
Long wavelengths are centrifugally stabilized if growth is 
slow compared to an epicyclic oscillation. Collapse re- 
quires Sff > i7, i.e. A < GS/ffi. The joint wavelength 
constraint /{GS) < A < GS/fP gives the Toomre in- 
stability criterion Qt ^ 1, that is satisfied throughout 
region III. 

The fastest growing modes in region III obey 



7III - O.Ty^/Qa , /Cfff = 0.5rs/Q^ 



(64) 



including finite thickness effects. Since 7111 ^ 1, collapse 
is faster than an orbital time scale. 

5.4. Across the Qt ~ 1 Transition 

The boundary between regions II and III has Qt ~ 1, 
which regains significance for Tg 3> 1 . The sharp increase 
in growth rates across the region II to III boundary (seen 
in Fig. [7]) is hardly surprising since standard GI stabilizes 
here. The presence of gas drag guarantees instability and 
gives complex behavior. Practical interest in this region 
is restricted to late stage, possibly second generation, 
planetesimal formation in gas depleted disks. 

Fig. [9] demonstrates that both the neutral mode (NM) 
and a density wave (DW) can collapse, competing for 
the fastest growth. For Qt — 1.05 only the NM grows 
and both DW solutions are damped at all wavelengths. 
For stronger self-gravity, DW growth appears at shorter 
wavelengths. By Qt = 0.85, the DW branch has faster 
growth. This competition generates wavelength discon- 
tinuities when tracking the fastest growth. In Fig. |31 the 
3 cm gas depleted {purple triple-dot dashed) curve shows 
such a switch near 20 AU. In Fig. [7] [left panel) this 
switching is evident in the winding /C = 1 curve (labelled 
"0" in log units). 

Fig. [To] demonstrates even more complex behavior 
when the self-gravity function J^vv turns negative, i.e. 



unstable to standard GI. Bifurcations appear where NM 
and DW modes suddenly change character. These in- 
triguing features reveal interesting dynamical behavior 
and suggest that care is required in numerical searches 
for fastest growing modes. 

5.5. Marginal Coupling 

We now consider the special case of Tg = 1, when or- 
bital and drag times are identical. This case is of special 
interest because drift times are fastest and many parti- 
cle concentra tion mechanisms are most efficient (CYIO, 
lYoudinl 120101 ). Furthermore numerical simulations are 
easier (not easy) without separate time scales. When 
Tg = 1 and Qa ^ 1 (i.e. the boundary between regions I 
and II) the limiting dispersion relation for the NM is 



7n 



Qt 



/C - (65) 



indicating stabilization by both effective pressure and dif- 
fusion as in Region II. The numerical coefficient, Cm — 
8/65, follows from our turbulence model [equation ([55)) ]. 
which should not be trusted to high precision due to com- 
plex dynamics in the Tg « 1 regime. The fastest growing 

modes obey 7inarg = /Cmaig/2 = CralQ\- 

Physical values for the growth of Tg = 1 modes are 



^grow,marg 



2.1 X 10-^ 



AU 



0.29 



a-4/T 



ring, mars 



lMa^4fTVR_ 



AU 



-i?Au , (66a) 
(66b) 
(66c) 



with a-4 = a/10~^. Despite marginal coupling, growth 
is still slower than dynamical and has ^ AU scales for 
standard disk parameters. 

We now extend the analysis to Qa — 1 and Tg = 1, 
the center of our parameter space. Solution of the full 
dispersion relation shows that density waves are damped 
and the fastest growing neutral mode has 7 w 0.1 and 
/C « 0.2. We make use of this result in comparing to 
simulations below. 

6. DISCUSSION 

6.1. Comparison to Simulations 

We now compare our linear instability with simula- 
tions that produce gravitation al collapse of Tg ~ 1 solids 
([Johansen et al.ll2007l l2009b[ ). We show that the non- 
linear gas dynamics included in the simulations gives 
faster growth than secular GI predicts. The relevant lin- 
ear growth rates for Tg = 1.0 are given in ii5.5l The fact 
that the simulations include solids with somewhat tighter 
coupling (smaller Tg) enhances the differences highlighted 
here. 

iJohansen et al.l (|2007l hereafter JetalOT) find grav- 
itational collapse occurs very rapidly, within a few 
orbits of introducing self-gravity. These simulations 
have turbulence generated by the MRI and by particle- 
gas interactions: both SI (streaming instabilities, 
lYoudin fc GoodmanI |2005( ) and vertical shear instabil- 
ities. Turbulent fluctuations in JetalOT have a Mach 
number T4ddy/cg ~ 0.05, which we characterize as a « 
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Figure 9. Growth rates (top panels) and wave frequencies (bottom panels), (in orbital units) vs. wavenumber IC for weak drag (ts = 30) 
and Toomre parameters: Qt = 1-05 (left panels) and Qt = 0.85 (right panels). The zero frequency neutral mode (NM, solid curves) is 
unstable at small K. Density waves (E)W, dashed curve) are damped for Qx = 1.05, but exhibit growth when Qt = 0.85 for K, near unity. 
Since J^yf > (grey dotted curves, bottom panels) all DW solutions would be stable in the absence of gas drag (Q-p = 0.85 is stable for 
standard GI due to a finite thickness correction). Analytic fits to the growth rates (dotted curves, top panel) from a Ts 3> 1 expansion work 
well except near J-\v ~ 0. 



0.05^ — 2 X 10^^. Their gas disk is reasonably massive 
with Qg ^ 15. The cohapse of sohds occurs not from a 
smooth background but from a ring-like clump generated 
by trapping in pressure maxima and Sl-clumping. Fig. 
2 of Jetal07 shows that this ring has a width of ~ 0.2hg. 
The mean particle surface density is roughly a factor of 
10 higher than the global average, giving Z ^ 0.1. This 
value is a generous assessment of the mean surface den- 
sity of the knotted ring. However, we now show that 
it is not generous enough to give agreement with linear 
theory. 

The linear growth rates presented above give collapse 



in a few orbital times when Qa ^ 1. This measure of 
self-gravity in turn requires Z ^ 0.5, from equation ([23]) 
and the values of a and Qg inferred above. This Z value 
is larger than the average value of the simulation clumps. 
Making the discrepancy worse, the linear theory requires 
an enhanced Z over an extended radial distance. The 
relevant dimensionless wavenumber /C ^ 0.2 corresponds 
to 



Ag 



lOvr. 



i0.2)hg 



1, 



(67) 



i.e. a wavelength that equals the gas scale height, wider 
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Figure 10. Similar to Fig. [9] but with Qt = 0.7. Stronger self- 
gravity gives J^vv < between 0.7 < < 2.0. Without drag, 
standard GI gives a growth rate itReal(\/— .^w) in this region {grey 
dashed curves, top panel). Gas drag affects the growth rates and 
introduce two bifurcations. At the first bifurcation (near K ~ 0.7), 
the unstable NM becomes an unstable DW, while the damped DW 
pair splits into a NM and a DW. At the second bifurcation (where 
K, ~ 1.2, the NM and unstable DW merge into an overstable DW 
pair, while the damped DW transitions to a damped NM. 



than the clumps in the simulations. Combining the mean 
surface density and width discrepancies gives an even 
more significant discrepancy in the total mass required 
for rapid collapse. 

This comparison shows that the collapse in simula- 
tions is faster and on smaller scales than in linear theory. 
A deeper understanding of how SI clumping interacts 
with self-gravity is needed. One factor is that SI (and 
perhaps other) c lumpin g can reduce turbulent diffusion. 
iJohansen et al.l ()2009bf ) show that turbulent stirring of 
particles is weaker than average in dense clumps. Less 
radial diffusion (lower effective a) will enhance secular 
GI and give better agreement with detailed simulations. 

6.2. Planetesimal Sizes and Solar System "Clans" 



A major goal of planetesimal formation models is to ex- 
plain the sizes of planetesimals. T his is useful as input for 
later stage m odels of coagulation ('Goldreich et al.| j2004l : 
iBromlev fc Kenvon 2006; Kenyon & Bromley 200^ and 
debris disks (Ken von fc BromlevI |2010j) , which in turn 
evaluate the viability of planetesimal formation models. 

Minor planets in the S olar System offe r by f ar the 
most detailed constraints. iMorbidelli et al.l ()2009() claim 
that the largest asteroids with diameters D > 100 km 
were "born big" by a collapse process. They argue that 
coagulation and coUisional evolution models alone can- 
not reproduce the size distribution of the largest as- 
teroids together with surface and shock age conditions. 
iWeidenschillinsJ ()2010D . however, argues that coagulation 
of small, _D ~ 0.1 km, planetesimals (with undetermined 
origins) could give better ag reement. 

In the outer Solar System, INesvornv et al.l ([2010') pro- 
pose that similar mass D > 50 km binary Kuiper Belt 
objects could form by late stage fission of a collapsing 
proto-pl anetesimal. The matc hing colors of binary com- 
panions (jBenecchi et al.|[2"009() supports this hypothesis. 

Collapsing planetesimals with equivalent Z? > 100 km 
form in the simulations described in H6.ll Since we do 
not model the non-axisymmetric fragmentation of secu- 
larly unstable rings, we do not make specific predictions 
of planetesimal sizes for comparison. However D 3> 1 
km is expected since the surface density when the ring 
fragments will greatly exceed the original S. 

We do predict that each secularly collapsing ring will 
produce a clar0 of fairly identical planetesimals. The 
solids in a ring will be well mixed during the slow initial 
collapse. Since a disk's chemical composition and mass 
density should vary smoothly with distance and time 
(and further will have sharp transitions at ice lines), each 
ring inherits unique conditions during its birth. Neigh- 
boring rings, or a subsequent generation, produce a plan- 
etesimal clan with different chemical composition and 
possibly sizes. 

The radial zonation of spe c trally distinct aster- 
oid classes (Gradic fc Tedescol Il982f l supports this 
fragmenting ring hypothesis, and is difficult to ex- 
plain with gradual coUisional evolution. The Kuiper 
Belt contains dynamical class es with different col- 
ors (iDoress oundira m et aTl I2008D and size distributions 
(|Bernstein et al. 20041) . These signs may merely point 
to dynamical mixing of a continuously variable forma- 
tion environment. More detailed study and evolution- 
ary modeling will help clarify if populations are distinct. 
This could strengthen evidence for clans of planetesimals 
from the same birth ring. 

6.3. Radial Drift and Metallicity Thresholds 

Radial drift is usually the dominant constraint on sec- 
ular GI growth times, but it might be over-restrictive. 
Growth over several drift times, idrift7 is possible 
because tdrift increases as the orbital radius shrinks 
[equation ([47]) ]. Moreover, sinc e drif t pilc-ups in- 
crease the surface densi ty of solids (jYoudin fc Shuil200^ 
lYoudin fc Chiand I2004D . collapse could accelerate over 
several tdrift- While we cannot include this sort of global 
evolution in our local analysis, we use drift pileups (and 

^ "Clan" purposely avoids "family," which already has a specific 
meaning. 
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other enrichment mechanisms, see iChiang fc Youdinl 
I2OIOD to motivate consideration of super-Solar Z values. 

Drift times would be much longer in special regions of 
the disk where the radial pressure gradient (and 77) de- 
creases. Pressu re ma xima halt radial drift and collect 
solids { W hipplill972| ). Enhancing Z in turn accelerates 
secular GI. Moreover the Keplerian gas motion in a pres- 
sure maximum has no headwind to stir particles. Thus 
pressure maxima are a potential goldmine for planetesi- 
mal formation. 

Unfortunately the existence and properties of pressure 
maxima remain uncertain. And appealing to special lo- 
cations must also explain why planetesimal formation ap- 
pears to have occured throughout the Solar System (and, 
with less detailed knowledge, throughout exoplanet and 
debris disk systems). A possible downside to pressure 
maxima is that they might be asso ciated with turbulence , 
as with MRI induced zona l flows (iJohansen et al.ll2009al ) 
or hydrodynamic vortices (|Lithwickll2009f) . Whether tur- 
bulent stirring offsets the advantages of pressure maxima 
is uncertain and may depend on particle size. 

We now consider turbulence driven by particle-gas 
interactions (for the standard case that radial pres- 
sure gradients do not vanish). Both vert ical shearing 
(IGoldreich fc WardlHOTl ICuzzi et aLlll993D and stream- 



ing ( Youdin fc GoodmanI |200'5[ lYoudin fc JohansenI 
[2007t IJohansen fc Y oudin' '200/) instabilities can drive 
midplane mixing. Wcidcnschilling (1995) famously ar- 
gued that vertical shear instabilities would prevent stan- 
dard GI. We show that secular GI is not impeded for 
sufficiently super-Solar Z . 

Particle-gas interactions stir the particle layer to a 
thickness h ^ r\R, a value supported by analytic theory 
CWeidenschiUing' 1995; Sckiva 1998; Youdin & Shu '2002) 
and simulation (Johanscii et al. 2009b; Lcc ct al. 2010). 
To achieve this thickness, turbulence adjusts to a level 

« T,(77i?//ig)2 « 5 X 10-7^^^ (68) 



using equation (|28|) . For secular GI to proceed we require 
from equation ([SS)) , which gives 
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assuming p > pg in equations (j32[) and (j45p . This thresh- 
old requires substantial enrichment, especially if gas de- 
pletion lowers F . It provides our only reason to favor 
secular GI in more massive disks. 

Howev er, particle dri v en tu rbulence may weaken at 
lower Z. lYoudin fc Sliiil (|200 2') predicted a threshold to 
suppress vertical stirring of Tg ^ 1 solids at 



Z > Zys 



r]R 
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(70) 



IJohansen et al.l (j2009b[ ) found a threshold for strong SI 
clumping of Z > 0.015 — Zys/3 with particle sizes 0.1 < 
Ts < 0.25 in 3D , vertic ally stratified simulations. 

iBai fc Stond (j2010a[ ) both confirmed the Sl-clumping 
threshold and showed that it goes away (or shifts to 
higher Z as they "only" w ent up to Z = 0.03) when 
smaller sizes are included. IBai fc Stoni (|2010bl ) show 
that the clumping threshold varies roughly linearly with 



the strength of radial pressure gradients, as in equa- 
tion dTOj). 

The existence of metallicity thresholds offers hope that 
Tg ^ 1 solids can overcome midplane stirring and un- 
dergo GI. More simulations are warranted, but the direct 
simulation of Tg ^ 1 solids is computationally intensive, 
especially in large domains and over the long time scales 
needed for secular GI. 

6.4. Validity of the Model Equations 

We now discuss the validity of our linearized equations 
of motion [eq. (ITT|) ] as a model of midplane particle dy- 
namics. 

Radial drift does not appear in (or arise from) our mo- 
mentum equations because a Galilean transformation re- 
moves steady drift and the locally constant radial pres- 
sure gradient. Basic disk quantities (S, f2, Tg) will 
change on the radial drift time scale, tdrift- We avoid 
this complication by conservatively restricting our atten- 
tion to modes that grow faster than tdrift- Drift speeds 
also cause a transition to non-linear, turbulent drag laws, 
addressed in appendix [K\ 

By using a single value of the stopping time, we do not 
consider particle size distributions. Since growth rates 
and wavelengths vary with Ts, it is worth considering 
whether a realis t ic dis persion of sizes could collapse co- 
herently. iWardI (|2000D demonstrated that a "bimodal" 
disk with two particle sizes is subject to secular GI. For 
a given Tg, a range of wavelengths around the fastest 
growing mode are unstable [see equation (|57|) and Fig. [3] . 
Thus a range of particle sizes can compromise and col- 
lapse at a common wavelength. Depending on the width 
of the size distribution, the growth rate will decline rel- 
ative to the mono-dispersive case. 

Using fluid equations to model t he dynamics of solids 
subjec t to gas drag was discussed bv lYoudin fc GoodmanI 
((2005h . The issue of the appropriate particle effective 
pressure and viscous forces to include for this type of 
particle fluid has not been rigorously explored. For this 
paper we use a standard relation between random veloc- 
ity and pressure [the second term on the RHS of equa- 
tion pib)) ]. Since mass diffusion is a stronger stabilizing 
influence, the form of the effective pressure should not 
be too significant. 

We neglect viscous forces, i.e. ter ms o f the form 
Dd'^u/djr and Dd'^v/dx^ in equations (|llb|) and (jllcp . 
respectively. It is po ssible for such viscous fo r ces to intro- 
duce instabilities, as lSchmit fc Tscharnuted (|1995D show 
for planetary rings. For our system however, drag forces 
dominate viscous forces. Thus we conclude that turbu- 
lence is more effective at diffusing solids than viscously 
transferring their momentum. Furthermore, the precise 
form of viscous forces is uncertain, so we drop them for 
simplicity. 

We address the neglect of detailed gas dynamics in the 
following section. 

6.5. Previous (& Future) Work on Secular GI 

IGoldreich fc WardI (|1973D noted that: "The frictional 
effect of gas drag does destabilize axisymmetric perturba- 
tions for wavelengths larger than [Ag]." Thus they were 
well aware of the existence of s ecular GI, but did not 
include it in their calculations. iWardI (| 19761 12000D de- 
rived and analyzed secular GI that is similar to this work, 
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but neglected the mass diffusion induced by gas turbu- 
lence. Indeed we ca n reproduce the dispersion relations 
in equation (A-22) oflWardI ()1976D and equivalently equa- 
tion (31) of iWardI (I2OOO0 bv setting Q^, = in our equa- 
tion (fT6| . lYoudinl (l2005ai rbl) analyzed a similar dispersion 
relation, adding a more sophisticated model for turbulent 
stirring, but still neglecting turbulent radial diffusion. 

These (and our) analyses model the dynamics of a sin- 
gle "fluid," the particles. Gas motions are in steady state. 
A two fluid analysis evolves the acrodynamically coupled 
motion of solids and gas. Spiegel (1972) considered the 
two fluid Jeans problem, i.e. self-gravity in a static, non- 
rotating, infinite medium. He showed that solids collapse 
relative to a static gas background on scales below the 
Jeans length. The Jeans length in a gravitationally sta- 
ble gas disk is quite large 



Aj 



0.9 



f3/4n9/8 



Fm 



1/4 



AU 



(71) 



supporting our assumption that long wavelength secular 
GI behaves in the single fluid limit. However the role of 
d isk dynamics remains unclear. 

iCoradini et al.l (jl981t ) considered two fluid GI in a thin 
disk. They too conclude that solids do collapse through 
a stationary gas over a broad range of wavelengths. A 
reanalysis of their model, including the effects of turbu- 
lent diffusion, would help clarify the upper wavelength 
cutoff to single fluid behavior and further explore two 
fluid phenomena. 

The two fluid GI model of iNoh et al.l ([l99l nu- 
merically evaluates the linear growth of global non- 
axisymmetric modes, whose behavior depends on reflec- 
tion at disk boundaries. The relevance of these modes 
for planetesimal formation is unclear at this time, but 
non-axisymmetric modes should not be forgotten. 

In the absence of self-gravity the two fluid model pro- 
duces streaming instabilities, i.e. SI (jYoudin fc GoodmanI 
|2005[) . Vertical speeds are r equired, thus height - 
integrated two fluid models like ICoradini et al.l (|1981D 
do not capture SI. The interaction of secular GI with 
the non-gravitational SI clumping would thus require a 
3D (or 2.5D for axisymmetric motion) two fluid model. 
The SI radial wavelength, Asi ~ TsTjR, does exceed the 
standard GI wavelength as 
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provided Tg is not very small. For the longer wavelengths 
of secular GI, the interaction with smaller scale SI clump- 
ing is less clear_; 

Goodman fc Pindoil (|2000( ) studied secular GI and 
pure drag instabilities in a (single) single fluid model. 
Their drag instability produces clumping and is analo- 
gous to SI (see CYIO). Drag instabilities can occur in 
their single fluid height-integrated model because they 
model drag as a surface density-dependent turbulent 
stress on the midplane layer. As their formulation differs 
from a stopping time based approach, a direct compari- 
son is difficult. Goodman & Pindor include viscous forces 
(that we ignore), offering evidence that they do not sup- 
p ress secular G I. 

iSekival (|1983D considers a different type of single fluid 
GI, that of a midplane layer with perfectly coupled gas 



and "dust" (ts ^ 1 solids). This model complements 
our treatment of collapse through the gas. Sekiya finds a 
buckling mode where the dusty gas incompressibly rises 
from the midplane when the total density exceeds '-^ pR 
[equation ([3l)]. The sedimentation of dust from this blob 
could then give planetesimals. However turbulent diffu- 
sion (which was not considered) should still be an ob- 
stacle to collapse. Thus the relevance of pn for small 
Tg remains unclear. A thorough investigation of the two 
fluid problem is needed to relate these different limiting 
cases. 

7. CONCLUSION 

We explore gravitational instability (GI) acting on the 
particle sublayer in protoplanetary disks. Gas drag in- 
troduces a class of secular (i.e. dissipative) GI that differs 
from standard (dissipation free) GI. We extend previous 
work by including mass diffusion due to turbulence. We 
focus mainly on the behavior of small solids, with stop- 
ping times Ts ^ 1, to study their possible gravitational 
collapse into planetesimals. We also consider Tg > 1 to 
understand the transition from secular to standard GI. 

We prove that secular GI always produces collapsing 
modes. This feature qualitatively differs from standard 
GI, which must satisly threshold conditions on Toomre's 
Qt ^ Ij or equivalently a critical, Roche-like density. 
While secular GI is formally ever-present, it can only 
produce planetesimals if collapse is more rapid than par- 
ticle radial drift and disk dispersal. The growth rates and 
wavelengths of secular GI depend sensitively on proper- 
ties of the gas and particle disks, which we briefly sum- 
marize. 

Gas drag is a double-edged sword. It introduces the 
possibility of secular GI, but gives progressively slower 
growth for smaller Tg ^ 1 by limiting infall to terminal 
speeds. Thus with other things (turbulence and Z) being 
equal, collapse is favored for larger solids and in the outer 
disk where gas densities are lower. 

Enhancing the disk "metallicity" Z promotes GI, both 
by increasing the strength of particle self-gravity and by 
slowing radial drift. A higher total (gas & particle) disk 
mass does not directly beneflt secular GI, because higher 
gas densities exert stronger drag (lower Tg). 

Turbulence slows growth, primarily by diffusively 
smoothing small scale collapse. Slow collapse occus over 
many orbits and thus takes the form of sheared-out rings 
that contract radially. These initially wide rings contain 
up to ~ O.IMq of solids, which will eventually fragment 
into lower mass planetesimals by a process not described 
here. 

Whether small solids can form planetesimals by secu- 
lar GI depends on what constitutes a "realistic" level of 
midplane turbulent diffusion, which we characterize as an 
a parameter. Our upper limits on a are low, ~ 10~^ — 
10^"^, depending on disk parameters, as shown in Fig. 21 

Our diffusive a can be much smaller than the more 
commonly used "viscous" (subscript added to dif- 
ferentiate) that measures accretion stresses and angu- 
lar momentum transport fjPringlei il981i) . Disk accre- 
tion may occur in active surface layers, with the mid- 
plane a quiescent dead zone where the MRI (magneto- 
rotational instability) does no t operate ()Gammi i Tl996l: 
iPerez-Becker fc Chiand 1201 If ). Furthermore, angular 
momentum transport could be dominated by magnetic 
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Maxwell stresses or waves that do not contribute to tur- 
bulent diffusion. This argument applies (e.g.) to waves 
launc hed in active layers that penetrate into the dead 
zone (jFleming fc Stonel[2nn8l : lOishi et al.l[200l . Finally 
the tendency of particles to concentrate in turbulence 
can counteract diffusion and reduce the effective a (see 
fjnm CYIO). We conclude that disk accretion, with ob- 
servationally constra ined values of 10"'' < < 0.1 
([Andrews et al.l[2010l ) , does not prevent planetesimal for- 
mation with significantly lower values of diffusive a. 

The strictest barrier to small particle GI is still 
posed by parti c le-driv en turbulence, as proposed by 
iWeidenschillingl (|1995f ). Overcoming this obstacle re- 
quires localized reductions in radial pressure gradients 
and/or enhancing Z to super-Solar values as discussed in 
^6.31 The strong correlatio n of ex oplanets with host star 
metallicity (jJohnson et aLl |2010D almost requires that 
planet formation have a strong Z dependence. 

Further progress can be made by better understand- 
ing the relation between secular GI and non-gravitational 
particle concentration by streaming instabilities and tur- 
bulent fluctuations. Dynamical theory and simulation 
complements interdisciplinary research on planet forma- 
tion that spans coUisional physics and experiments, small 
bodies in the Solar System, and astronomical studies of 
disks and exoplanets. 
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APPENDIX 

A. NONLINEAR DRAG LAWS 

Here we provide complete prescriptions to calculate the dimensionless stopping time Tg = fitstop, that gives a drag 
acceleration /^j-ag — — AV/istop on a particle of mass m — (47r/3)/9,a'^ moving at a relative velocity AV. We focus on 
turbulent drag forces which have a nonlinear dependence on relative velocity, as this case was not described in i|3.3.2l 
We also discuss the implications o f non-linear drag for linear stability analyses such as ours. 

For compact notation (following lAdachi et al.l[l976( ) we use the Knudsen number K = £g/a and the Mach number 
M = AV/cg. For the relative speed, we use the equilibrium drift speed between solids and gas. 



AV ^ r^QR "^ - , (Al) 

i.e. the vector sum of azimuthal headwind and radial drift, because it is much fast er than linear perturbation speeds. 
To include particle inertia , we could make the substitution (jNakagawa et al.lll986L see also equations [14] and [15] of 
lYoudin fc Goodmanl[200l 

rs^rs/(l + p/pg). (A2) 

This complication — making Tg depend on particle density, p, and thus turbulence and settling — is not justified for 
the current study. The correction vanishes for large Ts 3> 1 + p/pg, and thus has minimal effect (since large Ts is 
required in practice for turbulent drag). 
We assign Ts as follows: 

TEp if _ I < 



rsto = g^TEp if (I) M<K<^ 
-int ^ ^ (f )^/^sto if ^<X<(|)^M 



/^^2/5_ -r M ^ J. ^ i4\5 ■ (A3) 



Tturb = ^]g-TEp if K < 



Because the transitions themselves depend on via M, equation (jA3p is best solved numerically. We use integer 
coefficients to assure perfect matching at the transitions. The prescription for turbulent drag is consistent with the 
standard turbulent drag coefficient Cd ~ 0.44, since we can identify 

Cn^ - '""'^ - '\/.-^0-44 (A4) 

na^PgAV^ 3tstopPgAV inurhM 3-2003/5 ^ ^ 

Intermediate drag, given by Tint, describes the transition from viscous drag to fully developed turbulent wakes. 

We now consider corrections to the drag force in linear perturbation analyses when turbulent drag applies. In 
equation (jlip we used a drag acceleration of the form 

/drag = ^7 (A-5) 
''Stop 

where v = ux + vy is the linear perturbation velocity. While the total drag acceleration is always (anti-) parallel to 
the drift speed, for quadratic drag the perturbed drag force is not exactly (anti-) parallel to the perturbed velocity. 
Perturbed drag forces are twice as strong along the direction of the background drift. 
To modify equation ()A5P for quadratic drag, start with the total drag force 

f^,.,g = -j^\AV -v\iAV-v) (A6) 

^stop 

The stopping length £stop = I^V'l^stop is velocity-independent. The orientation of pressure gradient drift is 

and equation (IA2p would again give the inertial correction. The steady drag force is just — |AV|AV/4top- The linear 
perturbation is 



J dra 



stop 



2 sin2^\ , / . 2 " sin 2^ 



u cos u + V I X + f sm U + u I y 



(A8) 



We do not include this correction, as it would mainly apply to large Tg in the inner disk, not our primary focus. 
Generalizations to intermediate (not quadratic) drag and to include vertical motion are possible. 
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B. REGIONS IV & V: TOO DYNAMICALLY COLD 

We now consider two regions of our parameter space, which are unhkely to be of much astrophysical relevance. We 
include their description for completeness and to better understand the limitations of the mathematical model. 

B.l. Region V: Q'l'^ < < 1 

3/2 

Diffusion is exceedingly weak in region V. The condition ^ '^s corresponds to 

a « av ^ rlZ^IQl « 10-i«|^i?iu • (Bl) 

For comparison, this a is even lower than molecular viscosity, cknw ~ ^g/hg ~ 10~^^ R^^/F, except in the outer disk. 
In Fig. [21 the upturn in growth rates of the a = 10"-'^° curve beyond 30 AU is due to the transition to region V. 

In Region V free-fall gravitational collapse is countered my mass diffusion. Even though Tg ^ 1, collapse is so rapid 
that inertia exceeds drag forces. Balancing free-fall collapse [equation (l63l) ] with the diffusive time scale tdiff ^ ^^/Dp, 
gives the wavelength. Ay ~ D^/^ / {GSY^^ , and rate, sy ^ (GZ')^/'^/!)^/'^, or growth. An expansion of the dispersion 

relation confirms that 7y — l/Q^^ and K, — . 

With 7y ^ l/7"s, collapse is even faster than the (short in orbital terms) stopping time. This extremely rapid 
collapse introduces a physical inconsistency. The standard assumption that Tg <C 1 solids diffuse at the gaseous rate, 
D = Dg, likely breaks down. 

B.2. Region IV: Ql < 1/ts < 1 

Region IV is the large Tg partner of region V. The nominal requirement on a is again quite small. The more 
fundamental concern is that stirring by self- gravity would prevent a disk of Tg 3> 1 solids from being so quiescent. This 
was discussed in ^b.'6\ 

Putting these concerns aside, the behavior in region IV is identical to region V: gravitational collapse balances 
diffusion. This balance differs from region III, where the longer wavelength collapse is primarily opposed by pressure, 
not diffusion. The only change from the region V derivation is that D ^ Dg/r^ for Tg ^ 1. This gives 71V ~ {Ts/QaY^^ 
and /Civ ^ (t's/Qq)^^^- The collapse rate is much faster than dynamical with 7iy 3> Tg ^ 1, but this behavior does 
not appear in our disk models. 



